Transient and stationary behavior of the Olami-Feder-Christensen earthquake model 
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Using long-term computer simulations and mean-field like arguments, we investigate the transient 
time and the properties of the stationary state of the Olami-Feder-Christensen earthquake model as 
function of the coupling parameter a and the system size N. The most important findings are that 
the transient time diverges nonanalytically when a approaches zero, and that the avalanche-size 
distribution will not approach a power law with increasing system size. 
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I. INTRODUCTION 

The Olami-Feder-Christensen (OFC) earthquake 
model is probably the most studied nonconservative 
and supposedly self-organized critical (SOC) model. Sys- 
tems are called self-organized critical if they reach a sta- 
tionary state characterized by power laws without the 
need for fine-tuning an external parameter such as the 
temperature. Many researchers in the field agree on con- 
fining the term self-organized critical to those systems 
that are slowly driven and that display fast, avalanche- 
like dissipation events. This means that there is a sepa- 
ration of time scales, which can be interpreted as a way 
of tuning a parameter to a small value |2|. 

The prominent example for self-organized criticality is 
the sand-pile model by Bak, Tang and Wiesenfeld Q 
(BTW), where it can be shown analytically that the 
avalanche-size distribution is a power law, implying a 
scale invariance: Avalanches of all sizes are due to the 
same mechanism. The BTW model satisfies a local con- 
servation law, which can naturally lead to power laws 
0, , and without local particle conservation the model 
is not critical Q. The mechanisms leading to SOC in 
nonconservative systems are not yet well understood, and 
for the OFC model there is yet no agreement on whether 
it is critical at all. While some authors find critical be- 
havior when going to larger systems sizes and employing 
multiscaling methods 0,13, others interpret similar data 
as showing a breakdown of scaling |8( , and groups using 
branching-ratio techniques claim to find what they call 
almost criticality |g.llO|. 

Despite the simplicity of its dynamical rules, the OFC 
model shows a variety of interesting features that are un- 
known in equilibrium physics and appear to be crucial 
for generating the apparent critical (or almost critical) 
behavior. Among these features are a marginal synchro- 
nization of nei ghbo ring sites driven by the open bound- 
ary conditions and the violation of finite-size scal- 
es 0> 113 together with a qualitative difference between 
system-wide earthquakes and smaller earthquakes 
Also, small changes in the model rules (such as replacing 
open boundary conditions with periodic boundary condi- 
tions [l3|. introducing frozen noise in the local degree of 
dissipation [l4| or in the threshold values including 



lattice defects 0]), destroy the SOC behavior. Recently, 
it was found that the results of computer simulations are 
strongly affected by the computing precision |17| . and 
that the model exhibits sequences of foreshocks and af- 
tershocks 0,0. If energy input occurs in discrete steps 
instead of continually and if thresholds are random but 
not quenched, one finds quasiperiodicity combined with 
power laws |20j . The SOC behavior fully breaks down 
in OFC systems in one dimension [2l|, where only small 
and system-wide avalanches are observed. 

Since dynamics become extremely slow for large sys- 
tem sizes and for small values of the control parameter 
(implying strong dissipation) it is very difficult to obtain 
reliable results for the model based on computer simula- 
tions only. Thus, we find in the literature contradicting 
results concerning the transient time needed for the inva- 
sion of the 'self-organized region' from the boundary into 
the middle of the system, and concerning the avalanche- 
size distribution. While the transient time is found by 
some authors to scale with system size with an exponent 
depending on the level of dissipatio n ll ill . this exponent 
is found by others to be a constant [2q|. while still oth- 
ers find that above some critical degree of dissipation the 
invasion stops and never proceeds to the system's center 

m 

Similarly, the avalanche-size distribution is found ei- 
ther to be a power law with an universal exponent inde- 
pendent of the level of dissipation for large enough sys- 
tem sizes (however, different values for this exponent are 
reported in |'23| and in B, B), or a power law with a 
nonuniversal exponent |24| . Some authors found no 
power law at all above a critical degree of dissipation, but 
disagree on the value above which no power laws occur 
l2t| . Still other authors suggest that the dissipative 
OFC model is not critical at all (just like the random- 
neighbor version of the model |27| , which is a mean field 
approximation |28ll29ll3(]| ). but displays the new feature 
of being close to criticality, as mentioned above. If this 
is correct, only the conservative case leads to power laws 
in the distribution of avalanches [8|, |9| • 

By combining extensive computer simulations with an- 
alytical arguments, we will in this paper propose a phe- 
nomenological theory for the transient as well as the sta- 
tionary behavior particularly in the limit of large dis- 
sipation. The most important conclusions are that the 
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transient time diverges nonanalytically when the control 
parameter a approaches zero, and that the avalanche-size 
distribution will not approach a power law with increas- 
ing system size. 

The outline of the rest of this paper is as follows: In 
the next section, we present the definition of the model 
and explain the simulation algorithm. Then, we investi- 
gate the transient dynamics that brings the system from 
a random initial state to the stationary state as function 
of the system size and the model parameter. Section 
IIVI investigates the scaling behavior of the self-organized 
patches displayed by the system in the stationary state. 
The results flow into the interpretation of our simulation 
results for the avalanche-size distribution, which is stud- 
ied in Section [V] Finally, in Section IVII we summarize 
and discuss our findings. 



II. THE MODEL 

The OFC model originated by a simplification of the 
spring-block model by Burridge and Knopoff [3l| • To 
each site of a square lattice we assign a continous vari- 
able Zij £ [0, 1] that represents the local energy. Starting 
with a random initial configuration taken from a con- 
stant distribution, the value z of all sites is increased at 
a uniform rate until a site ij reaches the threshold value 
z t = 1. This site is then said to topple, which means 
that the site is reset to zero and an energy a x Zij is 
passed to every nearest neighbor. If this causes a neigh- 
bor to exceed the threshold, the neighbor topples also, 
and the avalanche continues until all zj-i < 1. Then the 
uniform increase resumes. The number of topplings de- 
fines the size s of an avalanche or 'earthquake'. The cou- 
pling parameter a can take on values in (0, 0.25). Smaller 
a means more dissipation, and a — 0.25 corresponds to 
the conservative case. Apart from the system size N, the 
edge length of the square lattice, a is the only param- 
eter of the model. Except for the initial condition, the 
model is deterministic. After a transient time, the sys- 
tem reaches an attractor of its dynamics. For periodic 
boundary conditions, the attractor is marginally stable 
and has a period of TV 2 topplings for all a 0, 01 • All 
avalanches have the size 1, and a site topples again only 
after all its nearest neighbors have toppled. Measured in 
units of energy input per site, the period is 1 — 4a. The 
behavior of the model is completely different for open 
boundary conditions, where sites at the boundary receive 
energy only from 3 or 2 neighbors and topple therefore 
on an average less often than sites in the interior. This 
leads to the formation of "patches" of sites with a sim- 
ilar energy, and this patch formation proceeds from the 
boundaries inwards. We are using open boundary condi- 
tions throughout this paper. 

Computer simulations of the model suffer from the long 
times needed to reach the stationary state for large N or 
small a. Most of the time is spent on searching for the 
site that will start the next avalanche, i.e. for the site 



with the largest value of z. Grassberger therefore used 
an algorithm that searches only among the sites with the 
largest values of z [l^. In our simulations, we used a 
different algorithm, based on a hierarchical search. The 
system size is chosen to be a power of 2. The system 
is divided into 4 boxes, each of which is again divided 
into 4 boxes, etc., down to the box on the lowest level, 
which consists of 4 lattice sites. Each box knows which of 
its 4 subboxes contains the site with the largest z value. 
Thus, the number of steps to find the site with the largest 
z value is log 2 N, since after an avalanche only those 
boxes have to be updated that have been affected by the 
avalanche. 



III. TRANSIENT TIME 

The transient time is the time needed for the patch for- 
mation to reach the center of the system. FigureHshows 
a system with N = 128 and a = 0.09 at three different 
times, the last snapshot being taken in the stationary 
state. 

One clearly distinguishes the patches close to the 
boundaries and the disordered inner part of the system, 
which behaves as if it was part of a periodic system. The 
time needed to establish a patchy boundary starting from 
a random initial configuration is very short, and virtually 
all of the transient time is needed to expand the patchy 
region to the entire system. The patches become larger 
with increasing distance from the boundary. 

The first ones to investigate the transient behavior 
were Middleton and Tang [TJ, who found that the tran- 
sient time increases as a power of N, with an exponent 
that depends on a. Later work on larger systems and 
for values of a larger than 0.15 by Lise |22j | found an 
exponent around 1.3, which does not depend on a. 

We will argue that the exponent does indeed depend 
on a, and that it diverges for a — ► 0. Figure [21 shows 
our simulation results for the transient time for different 
system sizes N as function of a. 

Each data point is based on the simulation of one sys- 
tem. Averaging over several initial conditions is not pos- 
sible because of the long computation times. The tran- 
sient time increases with increasing a and L, and it ap- 
pears to diverge for a — ► 0. For small a, one might 
therefore obtain the impression that the dynamics get 
completely stuck before the disordered block vanishes, as 
was suggested by Grassberger 0. However, we found 
no solid evidence and no good reason why this process 
should stop before the patches fill the entire system. The 
dash-dotted lines in Figure [3 are a fit with the Ansatz 
(|14f> . which is a generalized version of the result obtained 
in the following by using mean-field like arguments. 

We start with a local balance equation, which will lead 
us to an expression for the toppling profile as function of 
time. This consideration is similar to the one applied in 
|2l| to the one-dimensional model. Let iy be the mean 
number of topplings of site ij per unit time. If site ij top- 
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FIG. 1: Snapshots for a system with N = 128, a = 0.09 after 
10 5 , 5 x 10 5 and 8 x 10 5 topplings per site 



pies usually when zy is at the threshold (and not above) , 
tij must equal the mean amount of energy that this site 
obtains per unit time. For small values of a this assump- 
tion is well satisfied. Let g denote the rate of uniform 
energy input, and let a denote the average amount of 
energy passed to a neighbor during a toppling event. For 
small a, the value of a deviates very little from a, and 
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FIG. 2: Time measured in topplings per site until the inner 
block vanishes for N = 64 (circles), N = 128 (stars) and 
N = 256 (squares) as function of a; the lines correspond to 
the function T(a,N) = /(a)7V M(a) as derived in the text. 



in the later part of this calculation we will therefore re- 
place a with a. When discussing the size distribution 
of avalanches further below, we will see that with in- 
creasing system size the proportion of avalanches larger 
than 1 decreases towards zero, implying that the aver- 
age amount of energy passed to a neighbor approaches a 
even for larger values of a, and that most sites are ex- 
actly at the threshold when they topple, as was already 
observered numerically in |l7l l32| . The assumption that 
the value of a is constant throughout time and through- 
out the system is a mean-field assumption. Due to the 
approximations involved, we can expect that our theory 
makes predictions that are qualitatively correct, but that 
the quantitative features could be different. 
The balance equation reads 



U 



■ g + a{t i+ ij + ti-ij + ijj-i + t lJ+1 ) . (1) 



Now, we have to take into account the structure of the 
system during the transient time. The outer part consists 
of patches of different sizes [24| , and sites sitting in the 
same patch have to topple equally often for the patch to 
persist for a long time (which is observed by watching the 
system on the computer screen) . The value of tij depends 
therefore on the distance to the boundary, which can be 
measured in terms of the number of patches, x, between 
site ij and the boundary. (We ignore here the fact that 
the system has corners, which should not fundamentally 
change the argument. In any case, one could consider a 
system that is periodic in one dimension and open in the 
other, in order to avoid corners altogether.) Sites in the 
disordered block topple like in a system with periodic 
boundary conditions, i.e. they receive the same input 
from all four neighbors. For these sites we have therefore 
t = g + 4crf, or 



t = tn 



9 



1 - 4a 



(2) 
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In terms of the parameter x, the above balance equa- 
tion for the patchy part of the system becomes 



and approaches 1 (from above) for a 
condition 



0. From the 



t(x) =g + a(t(x - 1) + t(x + 1) + 2t(x)) , (3) 
or, in a continuum notation 

— tx-TT t(x)-£=0. (4) 
a dx z a 

The boundary conditions are t(0) = (x = signifying 
the non-existent neighbor of a boundary site) and t(d) — 
to, with d—1 denoting the index of the patch next to the 
disordered block. The solution of the balance equation is 
then 



t(x) = t 1 



sinh (n(d — x)) 
sinh (nd) 



(5) 



we obtain then 



d max (a,N) 



In 



N(Q-l) 
2qo 



lnQ 



2 go 
lnQ 



(9) 



(10) 



for large enough system sizes, qo is some constant (the 
extension of the patches of the first generation). The 
result for small a is therefore 



where k is given by k — y(I — 4a) /a. 

Next, we have to consider the advancement of the 
patchy structure into the inner part of the system. A 
site that is part of the inner block can become part of a 
patch only if the difference of its energy value z to that of 
its outer neighbor is less than a. This difference changes 
with time due to the different toppling rates. The patch 
next to the inner block topples less often than a neigh- 
bor of that patch, which is part of the inner block, the 
difference in the number of topplings per unit time be- 
ing to [sinh(rt)] / [sinh («d)], which is obtained from (J5J 
by inserting x = d — 1. The difference in the number of 
topplings per unit time is identical to the rate of change 
of the difference in the energy value z between the two 
neighbors. When this difference has increased by 1, it 
has taken any intermediate value (in steps of size a) and 
has therefore certainly assumed a value smaller than a. 
At that moment, the site of the inner block becomes part 
of the patch. The time (or number of topplings per site) 
needed to add an additional site to a patch is therefore 
proportional to 



n c (a, d) 



sinh nd 
sinh k 



In the limit of small a, n c (a,d) is given by 



n c (a, d) ~ exp 



(6) 



(7) 



which has to be summed over all patches, weighted with 
the mean size of each generation of patches. The total 
transient time is therefore 



T{a, N) 



ax(a,JV) 

E 

d=l 



l(d)n c (a, d) 



(8) 



with l{d) being the extension perpendicular to the bound- 
ary of a patch of type d. Below in Section IIVI we will 
see that 1(d) ~ Q(a) d ~ l , where Q is a function of a only 



with the exponent fi(a) = 1 + ^i~q7^ ■ Using the 
ansatz Q(a) — exp (/(a)), also motivated in Section Hvl 
with a leading term f(a) ~ Aa a and A and a positiv, 
yields 



and 



fj,(a) = 1 + 



N 



1 



Aa a+0 - 5 



T(a,N)~[—f(a)\ exp 



(12) 



(13) 



Inspired by this result of the mean-field theory, we ex- 
pect that the transient time is for small a given by an 
expression of the form 



T(a,N)~f(a)N 



fj,(a) 



(14) 



The data shown in Figure |21 agree with this expression. 
The numerical values of the parameter are A ~ 32.6 
and a ~ 1.262. f(a) was fitted in the form f(a) ~ 
exp [— V(a + 0.01)" + D], but any other expression could 
be equally valid. We would like to stress that this Ansatz 
was motivated by the mean field exponent fi for the lead- 
ing dependence on N, which we think mirrors the true 
behavior correctly. 



IV. CORRELATION FUNCTION AND 
CORRELATION LENGTH 

As has become clear from the previous section, the size 
distribution of patches as function of N and a and of their 
distance from the boundary is an important feature of the 
system. It affects not only the transient time, but also 
the avalanche size distribution, which will be discussed 
in the next section. 

We therefore investigate in this Section how the ex- 
tension of the patches in the directions parallel and per- 
pendicular to the boundary increases with the distance 
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from the boundary. For this purpose, we evaluate the 
correlation function 



(15) 



for a fixed distance i from the boundary for different 
times, starting again at a random initial configuration. 
We performed our simulations with systems that are pe- 
riodic in the direction of the second coordinate, i.e. site 
j + N is identical to site j. We chose N — 2 15 in order to 
obtain good statistics. The length of the system in the 
other direction was chosen just as large as needed, be- 
tween 48 (for small a, where the invasion front proceeds 
very slowly) and 512 (for large a). 

Figure |3 shows the correlation function for a = 0.08 
at distance 10 and distance 20 (measured in number of 
sites) from the boundary for three different times. One 
can see that at distance 10 the correlation function does 
not change any more with time, which means that the 
patch structure has been established at least up to this 
depth before the first measurement. We can furthermore 
conclude that the typical scale of the patches at a given 
distance from the boundary does not change any more 
when new patches are formed further inside. At dis- 





FIG. 3: Correlation function C(r) for a = 0.08 for a dis- 
tance d = 10 sites (top) and d — 20 sites (bottom) from the 
boundary at three different times (after 2570 (solid line), 5130 
(dashed line) and 12490 (dotted line) topplings per site). 

tance 20, we see that the correlation function builds up 



with time from zero to an exponentially decaying func- 
tion C(r) ~ e~ r /£. Figure U| shows the correlation length 
£ for a = 0.12 as function of the distance to the boundary 
for three different times. In the region where the patches 
are already present, £ increases as a power law in the 
distance from the boundary, and then falls down to zero. 
We see again that £ remains constant once the patches 
have emerged. The large fluctuations seen before the de- 
crease to zero occur in the region where the patches are 
just being formed. Due to large fluctuations in space, the 
averaging over the length of the system does not lead to 
a smooth curve for the system sizes used. 
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FIG. 4: Correlation length £ as function of the distance d to 
the boundary for a = 0.12 after 650, 1290 and 2570 topplings 
per site. 

Figure El shows the correlation length as function of 
depth for different values of a. The data are in good 
agreement with a linear increase of £ with the distance d 
from the boundary, but with a factor that decreases with 
decreasing a. However, we cannot rule out a power law 
£ ~ dP with an exponent r/ < 1 that increases with a. 

The linear (or power-law) increase of the correlation 
length together with the patchy structure leads to the 
following schematic picture (see Fig. 0) : 

The characteristic size of patches increases with dis- 
tance from the boundary. From one generation of patches 
to the next, the width and height of the patches increase 
with factors P(a) and Q(a) respectively. In the case 
r/ = 1, we have P = Q. (Of course, the patches at a 
given distance from the boundary do not all have exactly 
the same size, but a size of the indicated order of mag- 
nitude.) From snapshots of the systems, it is clear that 
P(a) and Q(a) increase with a. Furthermore, there must 
be a lower bound of 1 to both factors in the limit a — > 0. 
Thus, we can write Q(a) = exp(/(a)) with a monotoni- 
cally increasing function f(a) and f(0) — 0. The leading 
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FIG. 5: Correlation length £ as function of the distance to 
the boundary for a = 0.06,0.09,0.12 and 0.15 (from bottom 
to top) at the largest times simulated for a given value of a 
(solid lines); The dashed lines correspond to lines of slope 1 



P(a)p 




Q(a) qo 



FIG. 6: Schematic view of the system's structure: the width 
and height of patches increase with a power law in the distance 
to the boundary. Different generations of patches are coupled 
via n c , the increase in the size of the patches is P(a) parallel 
to the boundary and Q(a) perpendicular to it, starting with 
a size s = p <7o- 



dependence on a can be expected to be f(a) = Aa a with 
positive A and a. 

The correlation length in the ith generation of patches 
(counting from the boundary) is 



and the distance from the boundary is 

i 

3=1 



(16) 



(17) 



Based on this picture, we can write down an expres- 
sion for the size distribution of patches, which will be an 
important tool when discussing the size distribution of 
avalanches. A line at the distance d from the boundary 
cuts through ~ N/£h new patches of width £ and height 
fi £i/v through which the neighboring line does not 
cut, since there are ~ N sites along this line. The width 
distribution of patches is therefore given by 



N 



d(d). 



leading to 



N 

e 



(18) 



Transforming this into the size distribution np(s) with 
s ~ we obtain 



1 + 2TJ 

np(s) ~ Ns !+" 



(19) 



In the likely case that 77 = 1, we have an exponent —3/2 
in the size distribution of patches. Expressed in terms of 
P and Q instead of 77, the last equation becomes 



n P (s) ~ Ns 



(20) 



where P and Q depend on a. This expression can also 
be obtained directly from the recursion relation 



rip (sP(a)Q(a)) ds = / np(s) 



1 



P(a) 



ds , 



(21) 



where the integral is taken over one generation of patch 
sizes. 



V. AVALANCHE SIZE DISTRIBUTION 

Now we turn to the size distribution of avalanches in 
the stationary state. We made sure that the process of 
patch formation has reached the center of the system, 
before we evaluated the avalanche size distribution. 

In view of the results in Section IlIII we are now in the 
position to check how trustworthy the results reported in 
the literature are. As was already pointed out by Grass- 
berger |l2j . transient times are extremely long, and the 
first publications Q, S H E3 can have considered sta- 
tionary systems only for the largest values of a. 

It appears that many avalanche size distributions pre- 
sented in the last decade were actually obtained during 
the transient stage. We can check this only when the au- 
thors state how many initial avalanches they discarded 
for given N and a. Unfortunately, not all authors write 
how they decided if the system is in the stationary state. 
By observing statistical properties and comparing them 
at different times, one can be mislead to believe that the 
system has become stationary, although the advancement 
of the patches has only become very slow. Generally, the 
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larger a and the smaller the system size, the more likely 
is it that the published avalanche size distributions were 
obtained in the stationary state. For example, the re- 
sults published in 0,|2^| with L = 25, 45 were probably 
not taken in the stationary state for a below 0.2. Even 
Grassberger was evaluating avalanche size distributions 
during the transient stage in some parameter regimes. 

While taking small system sizes has the advantage of 
reaching the stationary state fast, they have the disad- 
vantage of being strongly affected by finite-size effects. 
It is therefore very difficult to predict the avalanche-size 
distributions in the thermodynamic limit. 

Figure shows avalanche size distributions for varying 
a with fixed N and for varying N with fixed a. The value 
of N in the top figure has been chosen small enough that 
the system could reach the stationary state even for the 
smallest value of a, which was 0.03. We can discern the 
following features: 

1. At least for value of a smaller than 0.17, the 
avalanche size distribution is no power law. A fit of 
the form n(s) ~ s -r(a)-a{a)\ns approximates the 
data much better than a pure power law. 

2. n(s) changes its shape with increasing N, implying 
that the system size affects the relative weight even 
of small avalanches, at least fo the system sizes con- 
sidered. This effect is stronger for smaller a. Only 
for the largest value of a is the main effect of the 
finite system size a rather sharp cutoff at iV 2 . 

3. The weight of avalanches of size 1 increases with 
increasing system size, while the weight of all larger 
avalanches decreases as 1 /N (see below) . 

In the following, we will explain these features based 
on the results obtained in the previous sections, and on 
what is known from literature. Described in words, the 
scenario is the following: Patches persist for a long time 
before they change their shape [241, due to an avalanche 
that enters the patch from outside [IJ, and patches fur- 
ther inside the system are rearranged less often. Large, 
patch- wide avalanches are mainly triggered at the bound- 
aries of the system. Whenever a patch-wide avalanche 
took place, there is a sequence of 'aftershocks' with de- 
creasing size according to Omori's law 0, and after a 
short time there occur mostly single topplings within a 
patch, until the next large avalanche comes from a patch 
of the previous generation. 

Let us quantify these statements. Analogous to the 
process of synchronizing neighboring sites discussed in 
Section ITTT1 neighboring patches also need a certain num- 
ber n c (a) of patch- wide avalanches in the patch closer to 
the boundary, before the inner patch experiences a patch- 
wide avalanche. This can be evaluated using Eq. (J5J for 
the situation that d — d max . We therefore obtain the 
recursion relation (compare l|21[l ) 



(P(a)Q(a)s) ds = / n pw (s) 



P(a)n c (a) 



ds . 



for the size distribution of patch- wide avalanches. If n c 
was independent of the generation number i, this would 
result in a power law n pw (s) oc A f s _T '°'" 1 with an expo- 
nent 



t{q) 



lnP(a)n c (a) 
lnP(a)Q(a) 



(23) 



For systems not too big, and for large enough a, there 
are only a few generations, and the approximation of a 
constant n c is not too bad. Evaluating Eq. JSJ for small 
a, we obtain the following result for n c that depends on 
the generation index i, 



n c (a, i) — exp( — =-) . 



(24) 



(see also equation |J7J|). Iterating Equation (|2*2*|l . we need 
to evaluate the product 



1 



n \p(a)n c (a,j)J \P(a) 



1 



exp 



2v^ 



which leads (using hxi ~ lns/lnPQ) to size distribution 
of patch- wide avalanches of the form 



where r(a) and a(a) are given by 



(26) 



r(a) = 
a(a) — 



ln(P(a)Q(o)) 
1 



LnP(a) - 



2y^ 



2^(ln(P(a)Q(a))) 2 



(27) 



(22) 



Now, we have to estimate the effect of 'aftershocks' 
on the size distribution of avalanches. These aftershocks 
will lead to an avalanche-size distribution that differs 
from that of the patch-wide avalanches. Aftershocks are 
avalanches that occur within a patch after a patch-wide 
avalanche. We assume that their size distribution is a 
power law with a cutoff at the size of the patch. This is 
motivated by the finding that systems that are dominated 
by one large patch display a power-law size distribution 
of avalanches (see also Therefore, we set 

n as (s\s') =s' T °s- r °8(s'-s), 

with n as (s\s') being the number of aftershock avalanches 
of size s in a patch of size s' , and with an exponent Th, 
which has a value around 1.8 (i.e., the value found in .7] 
for systems that have essentially one large patch). The 
size distribution of avalanches is then given by 



n(s) oc N J n pw (s')n as (s\s')ds 
= Ns~ To 

^ jy s -T(a)-cr(Q)lns 



/— r(a) — 1— a{a) In s' s 'to 



(28) 
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apart from a factor containing terms that depend on In s. 
Thus, the avalanche-size distribution is not a power law, 
but it has an exponent that depends logarithmically on 
s. As we have shown above, the data agree well with such 
a law. Figure |H| shows our results obtained for the coef- 
ficients a and r by fitting the avalanche-size distribution 
with the expression (|28fl . Although we can expect that 
the data are affected by finite-size effects particularly for 
small a, we see that the functions <j(a) and t{q) show a 
behavior that is in agreement with our expressions (|27ll : 
For small a, r decreases rapidly and will eventually be- 
come negative, while a tends to large positive values. 

The cutoff of the avalanche-size distribution is deter- 
mined by the size of the largest patch. As this size be- 
comes smaller with smaller a, the cutoff descreases also. 
Furthermore, since larger patches make a contribution 
to smaller avalanches via aftershocks, the effect of the fi- 
nite system size will be felt down to avalanche sizes much 
smaller than the largest patch. This is what is observed 
in the data. 

Finally, let us discuss the weight of avalanches of size 
1. After a patch- wide avalanche and the resulting after- 
shocks, a patch has single topplings (i.e., avalanches of 
size 1), just as a system with periodic boundary condi- 
tions, until a new patch-wide avalache comes from out- 
side. The total number of avalanches of size larger than 
1 per unit time is given by 

/>oo 

/ n(s)ds (x N . 

However, the total number of topplings per unit time is 
proportional to the number of sites in the system N 2 . 
We conclude, that only a proportion of the order 1 /N of 
avalanches has a size larger than 1. 

VI. CONCLUSION 

In this paper, we have shown that the transient time 
for the OFC model increases as function of the system 
size N and the coupling parameter a as T (a, N) ~ N a " , 
apart from corrections depending on a which do not af- 
fect this leading non-analytical behavior. This finding is 
in contrast to earlier predictions that the trasient time 
increases as a power law with system size, or that the 
transient time becomes infinite when a is smaller than 
some value. We obtained these results by performing a 
mean-field like calculation for the number of topplings 
per site and for the advancement of the patchy structure 
into the inner part of the system. 

Furthermore, by evaluating the correlation length of 
the energy values we found that the size of the patches in- 
creases as a power law with the distance from the bound- 
ary, leading to power law size distribution of the patches. 
Even if we assume that the size distribution of avalanches 



within a patch is a power law, we find based on the re- 
sults mentioned so far that the overall size distribution 
of avalanches is no power law, but has a logarithmic de- 
pendence in the exponent on the avalanche size, i.e., is 
of a log-Poisson form. This finding is supported by the 
simulation data and is valid at least for smaller a, where 
the system is not dominated by one large patch. 

We obtained our simulation results by using an effi- 
cient algorithm, however, the sharp increase of the tran- 
sient time with system size especially for small a made 
it impossible to study system sizes as large as necessary 
to see the true asymptotic behavior of the avalanche size 
distributions. 

Our findings are interesting for several reasons. First, 
the OFC model appears to show many features found in 
real earthquakes. As far as earthquake predictability j3|J 
or Omori's law [lsLIiII are concerned, this model appears 
to be closer to reality than others. If a is chosen above 
0.17, the avalanche size distribution agrees best with 
the Gutenberg-Richter law 0. Second, the OFC model 
demonstrates that apparent power laws need not reflect 
a true scale invariance of the system. We expect that this 
is true for many natural driven systems. Due to the dy- 
namics of the model, there occur avalanches of all sizes, 
however the mechanisms producing these avalanches are 
different on different scales. Large avalanches are mainly 
patch- wide avalanches, while smaller avalanches occur 
within patches during a series of foreshocks or after- 
shocks. Also, avalanches at different distance from the 
boundaries have different sizes. The observed "power 
laws" are thus dirty power laws, which appear like power 
laws over a wide range of parameters and over a few 
decades on the avalanche size axis, while the "true" an- 
alytical form is no power law. Third, the lack of a true 
scale invariance is accompanied by a decreasing weight 
of avalanches larger than 1 with increasing system size. 
This indicates again that the avalanche size distribution 
of the model does not approach some asymptotic shape 
with increasing system size, but that the weights of dif- 
ferent types of avalanches shift with the system size. This 
effect has most clearly been seen in one dimension, where 
the distributions split into a a dependent part at small 
avalanche sizes and a peak at sizes of order of the system 
size. Fourth, the extremely long transient times point 
to the possibility the some driven natural systems with 
avalanche-like dynamics are not in the stationary regime 
either. 
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FIG. 7: Size distribution of avalanches for different param- 
eter; top: system size N = 64 and a = 0.03,0.08,0.13,0.18, 
steeper curves correspond to smaller a; the s-axxis extends 
up to the total number of sites 4096; the distributions are 
normed on the total number of topplings; solid lines corre- 
spond to f(s) ~ s~ T ~ CTlns . middle and bottom: size distribu- 
tion for system sizes N = 64, 128, 256, 512; a = 0.09 (middle) 
and a = 0.17 (bottom); the distributions are divided by n(2). 
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FIG. 8: The coefficients a(a) (lower set of curves) and r(a) 
(upper set) as function of a as found by fitting the distribu- 
tions n(s) for N = 64, 128 and 256for those values of a where 
stationary systems were reached. 



